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We  are  interested  in  fitting  two-dimensional,  Gaussian 
conditional  Markov  random  field  (CMRF)  models  to  images. 

The  given  finite  image  is  assumed  to  be  represented  on  a 
finite  lattice  of  specific  structure,  obeying  a  CMRF  model 
driven  by  correlated  noise.  The  stochastic  model  is  char¬ 
acterized  by  a  set  of  unknown  parameters.  We  describe  two 
sets  of  experimental  results.  First,  by  assigning  values 
to  parameters  in  the  stationary  range,  two-dimensional  pat¬ 
terns  are  generated.  It  appears  that  quite  a  variety  of 
patterns  can  be  generated.  Next,  we  consider  the  problem 
of  estimating  the  unknown  parameters  of  a  given  model  for  an 
image,  and  suggest  a  consistent  estimation  scheme.  We  also 
implement  a  decision  rule  to  choose  an  appropriate  CMRF 
model  from  a  class  of  such  competing  models.  The  usefulness 
of  the  estimation  scheme  and  the  decision  rule  to  choose  an 
appropriate  model  is  illustrated  by  application  to  synthetic 
patterns.  Unilateral  approximations  to  CMRF  models  are 
also  discussed. 
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1.  Introduction 


Conditional  Markov  random  field  (CMRF)  models  have  many 
applications  in  image  processing  and  analysis;  for  instance, 
they  can  be  used  for  the  design  of  image  restoration  algo¬ 
rithms  [1-4] ,  for  characterization  of  textures  [4-7] ,  and  for 
image  coding  [4]  applications.  Apart  from  their  applications 
in  image  processing,  CMRF  models  are  an  active  research  topic 
in  the  statistical  literature  [8-9]. 

Typically,  an  image  is  represented  by  two-dimensional 
scalar  data,  gray  level  variations  defined  over  a  square  grid. 
One  of  the  important  characteristics  of  such  data  is  the  stat¬ 
istical  dependence  of  the  gray  levels  within  a  neighbor  set. 
For  example,  y(i,j),  the  scalar  gray  level  at  position  (i,j), 
might  be  statistically  dependent  on  the  gray  levels  over  a 
neighbor  set  that  includes  { (i-1 , j ) , (i+1 , j ) , (i , j-1) , (i , j+1) } . 
Due  to  the  lack  of  a  natural  preferred  direction  in  the  grid 
as  compared  with  the  discrete  time  interpretation  of  one-dim¬ 
ensional  Markov  processes,  the  classical  definitions  of  "past" 
and  "future"  are  not  generalizable  to  Markov  random  fields. 

Suppose  that  we  assume  that  y(i,j)  is  statistically  dep¬ 
endent  on  the  nearest  neighbors  along  the  east,  west,  north 
and  south  directions.  Then  the  Markov  property  in  the  cor¬ 
responding  CMRF  model  is  induced  by  requiring  that  the  con¬ 
ditional  probability  distribution  of  y(i,j),  given  the  values 
at  all  other  sites,  should  depend  only  upon  the  values  at  the 
four  sites  nearest  to  (i,j),  namely,  y(i-l, j) ,  y(i+l, j) , 
y(i,j-l)  and  y(i,j+l).  Wider  classes  of  CMRF  models  can  be 


considered  by  including  dependence  upon  the  values  at  more 
remote  sites.  For  instance,  corresponding  to  the  first  order 
Markov  model  mentioned  above,  we  can  have  a  second  order  CMRF 
model  that  includes  dependence  on  the  eight  nearest  neighbors, 
and  so  on.  In  this  paper,  we  are  primarily  interested  in 
Gaussian  Markov  fields,  i.e.,  the  conditional  distribution  of 
y(i,j)  given  all  other  lattice  observations  is  normally  dis¬ 
tributed  with  a  mean  that  is  a  linear  function  of  the  obser¬ 
vations  at  neighboring  sites  and  with  constant  variance 
(say)  v.  In  these  models,  the  observation  at  y(i,j)  is 
written  as  a  linear  weighted  sum  of  the  observations  over  the 
corresponding  neighbor  set  and  a  correlated  noise  sequence,  and 
is  characterized  by  a  set  of  coefficients  and  the  variance  of 
the  noise  driving  the  model. 

Prior  to  the  use  of  these  models,  two  problems  have  to 
be  tackled,  namely,  the  estimation  of  the  unknown  parameters 
and  the  choice  of  an  appropriate  neighbor  set  for  the  given 
image.  The  first  problem  has  received  some  attention  in  the 
literature.  Besag  [8]  developed  coding  schemes  for  the 
estimation  of  parameters  in  CMRF  models  and  maximum  like¬ 
lihood  (ml)  estimation  schemes  [9]  for  a  first  order  isotropic 
CMRF  model,  by  assuming  a  Gaussian  structure  for  the  noise. 

In  the  coding  method  the  grid  points  (or  sites)  are  divided 
into  two  sets  (say  A  and  B) ,  such  that  the  conditional  dis¬ 
tributions  of  gray  levels  associated  with  the  members  of  set  A 
given  the  observed  gray  levels  at  all  other  grid  points  are 


mutually  independent.  Once  the  forms  of  the  conditional  dis¬ 
tributions  are  specified  (except  for  a  few  unknown  parameters) , 
an  expression  for  the  conditional  likelihood  can  be  written, 
as  the  product  of  the  conditional  distributions  over  the 
members  in  set  A.  By  maximizing  this  likelihood  function, 
conditional  maximum  likelihood  estimates  of  the  unknown  para¬ 
meters  can  be  obtained.  One  of  the  main  disadvantages  of  this 
method  is  that  the  estimates  thus  obtained  are  not  efficient 
[9]  due  to  a  partial  utilization  of  the  data.  Though  by  com¬ 
plicated  coding  schemes  [10] ,  the  utilization  rates  may  be 
improved,  the  basic  problem  still  remains.  Also,  for  a  par¬ 
ticular  CMRF  model,  more  than  one  coding  scheme  can  be  realized, 
yielding  several  estimates  likely  to  be  highly  dependent.  Some 
arbitrary  methods  are  used  to  combine  these  estimates.  This 
coding  approach  has  been  used  for  the  analysis  of  geographical 
data  [8]  and,  very  recently,  in  texture  analysis  [7]  for  the 
case  of  binomial  variables. 

Unlike  the  cases  of  one-dimensional  time  series  models 
or  two-dimensional  unilateral  models,  deriving  an  expression 
for  the  likelihood  of  the  observations  poses  some  difficulties 
for  CMRF  models  that  include  dependence  on  neighbors  along 
all  directions.  This  is  due  to  the  fact  that  the  Jacobian 
of  the  transformation  matrix  from  the  noisy  variates  to  the 
observations  is  difficult  to  evaluate.  The  problem  of  eval¬ 
uating  this  determinant  can  be  avoided  by  making  assumptions 
about  the  representation  of  the  underlying  lattice.  Specif¬ 
ically,  by  assuming  representation  on  a  toroidal  lattice 


[4-6,9],  explicit  expressions  can  be  derived  for  the  de¬ 
terminant  term,  as  the  transformation  matrix  possesses  a 
block  circulant  structure  whose  eigenvalues  can  be  written 
down  explicitly.  We  use  this  representation  and  write 
explicit  expressions  for  the  likelihood  of  the  observations. 
Since  the  likelihood  function  is  nonquadratic  in  the  para¬ 
meters,  the  estimates  are  determined  by  using  numerical 
optimization  procedures  such  as  Newton-Raphson,  etc. 

An  alternative  estimation  scheme  that  yields  consistent 
estimates  for  the  parameters  is  given.  These  estimates  are 
similar  to  the  "least  square"  estimates  in  autoregressive 
time  series  models  (we  point  out  later  the  inappropriateness 
of  using  the  term  "least  square"  estimates  in  CMRF  models) . 
The  estimation  method  does  not  involve  computations  by  gra¬ 
dient  methods,  but  involves  inverting  an  mxm  symmetric  matrix 
for  a  CMRF  model  with  dependence  on  2m  neighbors. 

The  second  problem  considered  in  fitting  CMRF  models  is 
the  choice  of  appropriate  neighbors  in  images.  From  one 
dimensional  time  series  analysis  it  is  known  that  the  use  of 
an  appropriate  model  leads  to  good  results  in  forecasting 
and  similar  applications.  This  problem  has  been  considered 
in  the  literature  [4] ;  the  decision  rules  developed  are 
asymptotically  consistent,  transitive  and  parsimonious.  They 
are  based  on  the  corresponding  decision  rules  for  discrim¬ 
inating  between  different  autoregressive  models  [11] .  We 
implement  such  a  decision  rule  for  choosing  between  different 
neighbor  sets. 


The  usefulness  of  the  estimation  scheme  and  the  decision 
rule  for  the  choice  of  neighbors  is  demonstrated  by  application 
to  synthetic  patterns ,  the  underlying  true  model  of  thp  syn¬ 
thetic  patterns  being  known.  This  leads  us  to  the  problem 
of  generating  synthetic  patterns.  Computationally  elegant 
solutions  using  torus  representations  for  generating  synthetic 
patterns  obeying  CMRF  models  have  been  developed  in  [4] .  We 
use  this  approach  to  generate  a  number  of  two-dimensional 
Markov  patterns.  The  patterns  are  quite  varied  and  some  of 
them  possess  repetitive  periodic  patterns.  One  of  these  syn¬ 
thetic  patterns  is  used  for  illustrating  the  usefulness  of 
the  estimation  scheme. 

In  a  previous  paper  [12] ,  we  considered  fitting  another 
class  of  spatial  autoregressive  models,  the  so-called  sim¬ 
ultaneous  autoregressive  (SAR)  models.  Specifically,  we  sug¬ 
gested  an  approximate  ML  estimation  scheme  and  implemented  a 
decision  rule  for  the  choice  of  appropriate  neighbor  sets. 

The  CMRF  models  considered  in  this  paper  are  non-equivalent 
to  the  SAR  models  considered  in  [12] ,  in  that  given  a  SAR  model 
an  equivalent  CMRF  model  can  always  be  found  (usually  with 
more  parameters) ;  however,  the  converse  is  not  true.  For 
instance,  there  is  no  equivalent  SAR  model  for  the  first  order 
CMRF  model.  The  SAR  models  can  be  thought  of  as  higher  order 
Markov  models.  For  instance,  the  conditional  probability 
distribution  function  of  a  four-neighbor  SAR  model  with  dep¬ 
endence  on  the  east,  west,  north  and  south  neighbors  is  a 


function  of  the  eight  nearest  neighbors  and  the  second 
nearest  neighbors  on  the  east,  west,  north  and  south.  Some 
of  these  points  are  illustrated  later  using  synthetic  patterns. 

Another  topic  in  this  paper  is  unilateral  approximation  to 
CMRF  models.  Following  Besag  [8],  we  define  the  unilateral 
neighbor  set  of  any  site  (i,j)  in  the  positive  quadrant  to 
consist  of  those  sites  (k,£)  on  the  lattice  which  satisfy 
either  (i)ft<j  or  (ii)£=j  and  k<i .  Such  neighbor  sets  are  of 
interest  in  Kalman  filtering  of  noisy  images  [1,13].  We  take 
a  specific  case  of  a  first  order  CMRF  model  and  discuss  several 
unilateral  approximations.  The  patterns  corresponding  to 
the  first  order  CMRF  model  and  its  unilateral  approximations 
look  very  similar  in  structure. 

The  organization  of  the  paper  is  as  follows:  In  Section  2, 
we  consider  the  estimation  problem  and  suggest  a  consistent 
estimation  scheme.  The  implementation  of  the  decision  rule 
for  the  choice  of  appropriate  neighbors  is  discussed  in  Section 
3.  In  Section  4  experimental  results  are  given. 


2.  Estimation  of  Parameter  in  CMRF  Models 


Assume  that  the  observations  {y(s),s£ft}  have  zero 
mean  and  obey  the  CMRF  model  in  (2.1),  the  neighbor  set 
of  dependence  being  denoted  by  : 

y(s)  =  I  6,  .[y(s  +  (k,£))  +  y  (s  -  (k,Jt))]  + 

(k,  2.)  €Nf  k'£  (2.1) 

y'v  e(s)  ,  s£  n 

In  (2.1),  is  a  neighbor  set  such  that  if  (k,i)t  N^,  then 
( — k ,  —  2. )  ^  N.  ,  i.e.,  in  our  notation  =  {(0,1)  ,(1,0)}  means 
that  the  pixel  at  (i,j)  is  statistically  dependent  on 
neighbors  {  ( 0 , 1) ,  (0 , -1)  ,  (1 , 0) ,  (-1 , 0) } .  The  noise  sequence 
e(s)  is  correlated  with  zero  mean  and  unit  variance  and 
(0 ,  v)  are  unknown. 

Let 

q(s)  =  col[y(s  +  (k,£))  +  y  (  (s)  +  (-k,-2)  )  ,  (k,  2)fe 
Then  for  a  CMRF  model  we  have 

E  (e (s) q  (s) )  =  0  ,  Vs  (2.2) 

To  ensure  homogeneity,  the  coefficients  {0^  ^ ,  (k ,  2)jtN-^ } ,  must 
obey 

U ' i j  (0)  >  0  (2.3) 

where 

u’  •  .  (0)  =  (1-29T4>.  .)  (2.4) 

and 

4) .  .  =  col  [cos  2tt  (  ( i—  1 )  k  +  ( j-i)  2)  ,  (k,  Jl)fcN,  ]  (2.5) 

3  M 

For  a  finite  image,  and  a  CMRF  in  which  neighbors  in 
all  direction  are  considered,  some  of  the  neighbors  for  the 


boundary  sites  are  not  defined.  Several  different  assump¬ 
tions  can  be  made  regarding  the  distribution  of  these  pixels 
[4]  [14].  In  this  paper  we  assume  that  the  image  is  folded 

on  a  toroidal  lattice  such  that,  for  all  (i,j)tft, 

y[(s)  +  (i1,j1)]  =  y[(s)  +  (i.^-1 , j  1~1) mod  M  +  1]  (2.6) 


This  assumption  ensures  that  all  the  relevant  neighbors  of 
any  y(s)  belonging  to  the  finite  image  are  well  defined. 
Letting  e  be  the  lexicographic  ordered  array  of  {e(s),s  (2), 
(2.1)  can  be  rewritten  as 

H  (0)  y  =  /v  e  (2.7) 

where 


0  =  col  [8k  t,  (k,  NjJ 
H  ( 6 )  is  a  block  circulant  matrix 


Hl,l  H1 , 2 

•  •  • 

hi,m 

H  (9 ) 

= 

H1,M  Hl,l 

Hf , 2  '  •  * 

hi,m-i 

•  •  • 

_K1, 2 

Hi , 1  _ 

and  each 

component  matrix  H 

.  is  circulant. 

r  J 

For  i 

when 

=  (  (0 

,  1 ) ,  ( 1 , 0 ) }  we 

have 

Hl,i 

=  circulant  (1,8^ 

,°,' • • ,0O,-1) 

(2.8) 


^  2  =  circulant  (0^  g,0,...,0) 


M  =  circulant  (9_^  q,0,...,0) 


and  H.  . 

r  J 


0,  j  1,2, .  .  .  ,M 


Given  an  image,  we  are  interested  in  fitting  a 


CMRF  model  to  it.  This  problem  has  received  some  attention 
in  the  literature.  In  [8],  Besag  developed  coding  schemes 
for  estimation,  and  ML  schemes  were  used  in  [9]  for  simple 
CMRF  models.  As  mentioned  before,  the  coding  approach 
yields  inefficient  estimates  due  to  a  partial  utilization 
of  data.  ML  estimates  can  be  derived  by  assuming  some 
appropriate  distributions  for  the  noise  sequence  {e(s),s£Q}. 
Assume,  for  instance,  that  e(s)  is  distributed  normally  with 
zero  mean  and  unit  variance.  Using  this  assumption,  the 
log  likelihood  function  in  p(y)0,v)  can  be  written  as 

in  p  (y  |  0  ,  v)  =  in  det  H(0)  -  (M2/2)  in2irv  -  -  yTH(0)y  (2.9) 

2v~ 


From  the  theory  of  block  circulant  matrices, 

in  det  H (0 )  —  Z  Any*  . (0)  (2.10) 

(i,j)6ft  D  ~ 

Using  (2.4),  (2.9)  and  (2.10)  we  have 

in  p  (y  !  0  ,  v)  =  Z  in  (1—2  .  .)  -  (M2/2)in  2itv 
~  ~  (i,j)€ft  ~~13 

Since  H(0)  is  a  block-circulant  matrix,  it  is  diagonalized 
by  a  two-dimensional  FFT,  with  eigenvalues  u'^j(9).  Using 
this  fact,  we  have 


yTH(9)y  =  Z  ||  z  (A  .  . )  ||  (1-26T<K  .) 
~  ~  o  ~ 


where 


z(Aij> 


are  the  Fourier  transforms  of  y(s). 


sfefi  and 


X..  =  ( 2n  ( i-1)  ,  2n  ( -j-i )  )  .  Substitution  of  (2.12)  into  (2.11) 
13  M  M 

gives 

£n  p(y  1 0 ,  v)  =  Z  In  (1-20T$..)  -  (M2/2)  £n  2ttv 

'■  ~  o  ~  ~1J 

(2.13) 

_1  E  ||z  (A  )  |[2  (1-26T4>  ) 

2v  Q  13  ~  '-13 

Note  that  the  contribution  of  the  exponent  term  of  the 
likelihood  function  is  linear  in  0,  unlike  the  case  of  SAR 
models  [12],  where  this  term  is  quadratic.  Consequently, 
the  notion  of  least  squares  estimates  does  not  extend  to 
CMRF  models.  Numerical  optimization  procedures  such  as 
Newton-Raphson  can  be  applied  to  obtain  ML  estimates.  The 
ML  estimates  thus  obtained  are  asymptotically  consistent 
and  efficient.  However,  the  computation  of  ML  estimates 
might  become  expensive  due  to  the  calculation  of  the 
gradients  of  the  log  likelihood  function.  An  approximate 
ML  estimation  scheme  for  CMRF  models,  without  substantial 
sacrifice  in  accuracy,  is  considered  elsewhere  [15] .  To 
avoid  these  computational  difficulties,  we  suggest  the 
following  estimation  scheme:  consider  the  estimates 

6Q  =  [Z  q(s)qTs)]  (E  q(s)y(s))  (2  14) 


and  :  _  1  E  (y(s)  -  0^  z (s) ) 2  (2.15) 

vo  ~  m2  n  ~ 

A 

The  easily  computable  estimate  fig  is  asymptotically  consistent, 
as  we  will  now  prove. 


Theorem  1 :  Let  y(s),sfcft  be  the  set  of  observations  obeying 


the  CMRF  model  (2.1).  Then  the  estimate  6^  is  consistent. 
Proof :  From  (2.1)  and  (2.14), 

e  =  [Z  q(s)qT(s)]_1  (Z  q(s)(0Tq(s)  +  /ve(s)) 

~u  ft  ~  ~  ft  ~ 

=  0  +  [Zq(s)qT(s)]‘1/ve(s)  (2.16) 

ft  " 


Equivalently, 


[Z  q(s)qi(s)]  (9.  -  6)  =  /ve  (s) 

—  -N.  --  ~  ~ 


(2.17) 


Since  the  matrix  Z  qfs)ql(s)  is  positive  definite,  and 

a  ~ 

E(q(s)e(s))  =  0  with  probability  1,  the  assertion  of  the 
theorem  follows. 

A 

However,  the  estimate  0q  is  not  efficient  since  it  is 
not  equal  to  the  ML  estimate  obtained  by  maximizing  (2.13). 
In  this  paper  we  use  the  estimate  in  (2.14)  and  (2.15)  for 
simulation  studies. 


3.  Model  Selection 


We  briefly  discuss  the  decision  rules  for  the  choice 
of  appropriate  CMRF  models.  The  importance  of  appropriate 
model  selection  was  illustrated  in  the  case  of  SAR  models 
[12] ,  [15] ,  using  synthetic  patterns.  Specifically,  it 
was  shown  that  the  quality  of  reconstruction  varied  con¬ 
siderably  depending  upon  the  underlying  model.  We  can  ex¬ 
pect  a  similar  situation  with  respect  to  CMRF  models.  The 
problem  of  choosing  an  appropriate  CMRF  model  becomes 
difficult  due  to  the  rich  variety  of  possible  model 
structures.  Suppose  we  have  an  original  image,  say  a  64  x 
64  window  from  one  of  the  Brodatz  textures,  and  we  fit 
different  CMRF  models  to  the  texture.  It  can  be  argued 
that  by  visual  inspection  of  the  reconstructed  patterns 
corresponding  to  the  different  fitted  models,  a  decision 
can  be  made  regarding  the  appropriate  model.  However,  there 
are  several  objections  to  this  procedure.  First,  the 
decision  rule  is  subjective  and  no  quantitative  measure 
of  possible  error  in  the  decision  is  given.  More  signifi¬ 
cantly,  the  reconstructed  patterns  corresponding  to  an 
original  RF  model  and  another  RF  model  which  includes  the 
original  model  and  some  extra  neighbors  look  very  similar. 
Hence,  a  decision  based  on  visual  inspection  is  unreliable. 
Also,  given  an  arbitrary  image  pattern,  no  original  to 


compare  with  being  available,  it  should  be  possible  to 
choose  on  a  quantitative  basis  an  appropriate  model  from 
a  family  of  such  models.  In  the  context  of  this  paper, 
different  CMRF  models  should  be  interpreted  as  representing 
different  neighbor  sets  N^. 

The  problem  of  choosing  appropriate  CMRF  models  has 
received  some  attention  in  the  literature.  The  possible 
approaches  are  using  pairwise  hypothesis  testing  procedures 
[7-8],  Akaike's  information  criterion  (AIC)  [16],  and  the 
Bayes  approach  [11]  [17].  The  main  criticisms  of  the 

pairwise  hypothesis  testing  approach  for  CMRF  model  selection 
are  that  the  resulting  decision  rules  are  not  transitive 
[11],  and  the  decision  rules  are  not  consistent,  i.e., 
the  probability  of  choosing  an  incorrect  model  does  not 
go  to  zero  even  as  the  number  of  observations  goes  to 
infinity.  When  coding  schemes  are  used,  more  than  one 
coding  scheme  results  for  the  same  data.  It  might  well 
turn  out  that  the  null  hypothesis  may  be  rejected  on  some 
coding  schemes  but  accepted  on  some  other  coding  schemes 
for  the  same  data;  some  ad  hoc  methods  leading  to  very 
conservative  decisions  [7-8]  are  used  to  overcome  this 
problem.  When  it  is  required  to  choose  between  a  first  order 
CMRF  model  and  a  second  order  CMRF  model,  to  ensure  that  the 
likelihoods  are  comparable,  the  coding  scheme  corresponding  to 
the  second  order  CMRF  model  should  be  used  for  the  first  order 
CMRF  model  also,  causing  further  loss  of  efficiency  in  est¬ 
imation  for  the  first  order  model.  Also,  frem  a  philosophical 
point  of  view,  the  pairwise  hypothesis  testing  procedure  is 


not  suitable  for  choosing  from  a  family  of  different  CMRF 
models.  This  methodology  is  well  suited  for  classical  prob¬ 
lems  such  as  drug  testing  where  the  distribution  under  the 
null  hypothesis  that  the  drug  is  not  effective  is  specified 
and  the  consequences  of  rejecting  the  null  hypothesis 
when  it  is  not  true  (type  I  error)  are  crucial.  On  the 
other  hand,  in  the  model  selection  problem,  we  are  only 
interested  in  choosing  a  model  that  can  reasonably  account 
for  the  observed  statistical  phenomena.  When  the  problem 
is  posed  in  such  a  way  that  the  null  hypothesis  is  a 
specified  CMRF  model,  and  the  alternative  is  its  negation, 
the  consequences  of  type  I  error  are  not  serious.  In 
other  words,  we  cannot  take  a  subjective  view  of  any 
specific  hypothesis,  or  equivalently  any  specific  CMRF 
model,  as  is  implicitly  required  in  pairwise  hypothesis 
testing  procedures. 

The  model  selection  problem  comes  under  the  category 
of  multiple  decision  problems.  A  method  that  is  well  suited 
for  this  problem  would  be  to  compute  a  test  statistic 
for  different  models  and  choose  the  one  corresponding  to 
the  minimum.  The  AIC  criterion  and  the  Bayes  method  are 
two  such  procedures.  The  AIC  statistics  can  be  formulated 
from  the  expression  of  the  log  likelihood  function  given 
in  Section  2;  the  best  model  is  the  one  which  minimizes 
the  AIC  statistic.  The  method  gives  transitive  decision 
rules  but  is  not  consistent  even  for  one-dimensional 


autoreqressive  models  [18] .  Hence  it  is  not  advisable 
to  use  the  AIC  method  for  CMRF  model  selection. 


We  formulate  the  problem  and  suggest  a  decision 
statistic  that  is  consistent,  transitive  and  parsimonious. 
The  actual  derivation  can  be  done  by  using  standard  Bayes 
decision  theory,  as  was  done  for  SAR  models  in  [17] . 
Suppose  we  have  three  sets  N}1'N12'^13  °f  neighbors  con¬ 
taining  m1 members,  respectively.  Corresponding 
to  each  N.  ,  we  write  the  CMRF  model  as 

iq 


y  (s)  =  L 

(k,  fc)€N. 


Og  k  e  [y(s  +  (k,Z)  +  y (s- (k, Z) ]  + 


v^Tc  (s)  , 

°q,k,e+  °'  (k'^€Nlq'  'Vq>0'  =  1,2,3 

Then  the  decision  rule  [4]  for  the  choice  of  appropriate 
neighbors  is:  choose  the  neighbor  set  N.  *  if 


q*  =  arg  min  {C  } 

_  q  q 


(3.2) 


where  *  '  2  -> 

C  =  -25  s?  n  (1-0  ({>..)+  M  inv  +  2m  £n  (M  )  (3.3) 

q  n  -q~qii  q  q 


where 


!q  *  C0l[nq,k,l  ' 


.  .  =  Col  [cos  2 ir  (  ( i - 1 )  k  +  ( j-1)  £)  ,  (k,  £)eM.  ] 

-  4  •  1  r  j  M  -^4 


The  model  selection  procedure  consists  of  computing 
for  different  models,  and  choosing  the  one  corresponding 
to  the  lowest  C^.  The  factor  2  appears  in  the  third 
term  of  (3.3)  to  account  for  the  fact  that  in  our 
notation  m  refers  to  the  number  of  symmetric  neighbors. 


4.  Experimental  Results 


We  describe  the  results  of  some  experiments  regarding 
the  generation  of  two-dimensional  patterns  obeying  known  CMRF 
models  and  estimation  schemes  developed  in  Section  2. 
Experiment  1:  Synthetic  generation  of  two-dimensional 

Markov  patterns. 

From  (2.7)  we  have 

H  (0 )  y  =  /ve  (4.1) 


where  e  is  a  zero  mean  correlated  noise  sequence  with  correla- 
T 

tion  matrix  E  (ee  )  =  H(0).  The  synthetic  generation  is  then 
done  by  assigning  some  arbitrary  values  in  the  stationary 
region  to  0,  and  using  the  correlated  noise  vector  e.  The 
sequence  (e(s)}  may  be  generated  by  using  DFT  computations 
[19],  and  y  can  be  generated  by  inverting  the  matrix  H(0). 
Since  H(0)  is  a  block-circulant  matrix,  Fourier  computations 
can  be  used  for  generating  y.  An  alternative  procedure 
developed  in  [4]  is  used  in  this  paper.  Before  proceeding 

2 

further,  we  define  the  following  quantities:  denote  the  M 


Fourier  vectors  f  ^  ,  i  <  i,j*M  by 

f.^  =  col  [tj  ,  XAtj  ,  . . . ,  X^rltj]  , 

tj  =  col  [1  ,X  j  , A  1]  ,M- vector 


A.  =  exp  [  2tt  (i— 1)  ] 
M 


The  synthetic  generation  scheme  then  is  as  follows: 


y  =  Z  (f .  .x.  .//i?.  .  (0)  )  +  ctl 

1  o  -*1]  xy  ~ 


(4.2) 


where 


x .  .  =  /v  E  f .  .  w 
13  2  9  3  ~ 


M 

1  =  (1,1,... ,1)  ,  M^-vector 

and  u  ’  i j  ( 0 ) #  lsi,jim,  are  defined  in  (2.4)  and  w  is  an  identical 
and  independently  distributed  noise  sequence  of  known  distribu¬ 
tion.  For  a  proof  that  such  a  y  indeed  obeys  underlying  CMRF 
model ,  see  [4 ] . 

We  generate  the  vector  w  using  a  Gaussian  pseudo-random 
number  generator,  generate  its  Fourier  sequence  {x^}  by  a  two-* 
dimensional  FFT ,  and  finally  use  (4.2).  Sixteen  such  64x64 
images  were  generated  using  CMRF  models  with  different  neighbor 
sets  and  parameters.  The  grayscale  values  of  the  images  were 
scaled  to  lie  in  the  range  0-63.  Alternatively,  by  multiplying 
the  value  of  v  used  by  an  appropriate  constant  the  same  patterns 
are  obtained  without  scaling.  The  details  of  the  models  are 
given  in  Table  I  and  the  corresponding  images  are  shown  in 
Fig.  1.  It  can  be  seen  that  the  patterns  generated  are  quite 
varied  and  some  of  them  look  similar  to  real  textures.  Contrary 
to  the  existing  belief  [20]  that  spatial  autoregressive  models 
are  incapable  of  exhibiting  the  local  pattern  replication 
attribute  considered  an  essential  ingredient  of  texture,  some 
of  the  windows  do  exhibit  repetitive  patterns. 


We  use  matrix  notation  in  referring  to  windows  of 
images.  Diagonal  neighbor  sets  seem  to  induce  diagonal 
patterns,  as  in  the  windows  in  positions  ( 1 , 2) ,  ( 1 , 3)  , ( 2 , 2) 
and  (2,3).  The  pattern  in  the  window  (2,1)  corresponding 
to  a  CMRF  model  with  =  {  (0 , 1) ,  (1 , 0) ,  (-1 , 1) ,  (1 , 1) ,  (0 , 2) , 

(2,0) }  can  also  be  generated  by  a  SAR  model  as  shown  in 
[12]  with  N  =  {  (0,1)  ,  (0,-1)  ,  (-1,0) ,  (1,0)  }  and  e_1  0  =  ®1  0  = 
-.12,  0g  _i  =  6q  i  =  .28  using  the  procedure  outlined  in  [4-5]. 
Note  that  this  pattern  and  the  one  in  window  (2,1)  of  Fig.  1 
in  [12]  are  very  similar.  Likewise  the  diagonal  pattern  (2,2) 
generated  by  a  CMRF  model  with  =  {  (-1 , 1) ,  ( 1 , 1 )  ,  ( 2 , 0) ,  (0 , 2 ) , 
(2,-2) ,  (2, 2)  )  can  also  be  generated  by  a  SAR  model  with  N  = 

{  (-1,1)  ,  (1,-1)  ,  (1,1)  ,  (-1,-1)}  ,  Q_ltl  =  B1'_1  =  -.14  and 
0^  ^  =  .28.  However,  the  patterns  in  windows  (1,1) 

and  (4,4)  corresponding  to  a  CMRF  model  with  =  {  (0 , 1)  ,  (1 , 0)  } 

cannot  be  generated  by  any  finite  parameter  SAR  model. 

The  role  played  by  adding  nearest  neighbors  can  be 
illustrated  using  patterns  (4,2),  (4,3)  and  (2,1).  Window  (4,2) 
corresponds  to  =  {  ( 1 , 0) ,  ( 0 , 1) ,  (2 , 0)  ,  ( 0 , 2)  }  and  produces 

horizontally  oriented,  macrostructured  strip  patterns.  By 
adding  the  symmetric  neighbor  (1,1)  a  diffused  version  in 
window  (4,3)  is  produced.  When  an  extra  symmetric  nearest 
neighbor  (1,-1)  is  added  we  obtain  a  more  microstructured 
pattern  resembling  water.  By  adding  similar  neighbors  to  the 
model  in  (4,1),  vertically  oriented  patterns  can  be  produced. 


The  role  of  parameter  values  in  the 
structure  of  patterns. 


Experiment  2 : 

To  illustrate  the  role  played  by  the  coefficients 
in  generating  the  two-dimensional  patterns,  we  consider  the 
pattern  corresponding  to  the  CMRF  model  =  {  (-1,1) ,  (1,1)  } , 
a  =  30-000,  v  =  1.1111,  0_^  ^  —  .28  and  0^  ^  =  -.14.  The 
values  of  the  parameters  tried  are  given  in  Table  II  and  the 
resulting  pictures  in  Fig.  2.  Note  that  as  the  parameters 
are  varied,  the  basic  pattern  is  still  retained  but  the 
"busyness"  of  the  pattern  is  varied. 

All  the  patterns  considered  thus  far  were  generated 
using  the  same  pseudorandom  number  generator.  As  illustrated 
in  Figs.  4  and  6  of  [12]  for  SAR  models,  different  psuedorandom 
sequences  produce  very  small  perturbations  in  the  patterns. 

Experiment  3: 

To  test  the  usefulness  of  the  estimation  scheme 
and  the  choice  of  appropriate  neighbor  sets,  experiments  were 
done  with  one  synthetic  pattern.  The  true  model  used  to 
generate  the  test  pattern  corresponds  to  =  {  (-1 , 1) ,  (1 , 1)  } , 
a  =  30,00,  v  =  1.1111,  0_^  ^  =  -.14,  0^  ^  =  .28.  Using  this 
CMRF  model,  the  synthetic  image  (1,1)  in  Fig.  3  was  generated. 

But  for  correct  inference  purposes  regarding  the  estimation 
schemes,  the  original  window  (values  not  scaled  for  display 
purposes)  was  used.  For  estimation  of  the  parameters,  the 

i 

f 


sample  mean  of  the  window  was  subtracted  and  (2.14)  and  (2.15) 
were  used.  The  test  statistic  c  in  (3.3)  was  also  computed. 

The  actual  values  of  the  estimates  corresponding  to  different 
neighbor  sets  are  given  in  Table  III  and  the  corresponding 
reconstructed  images  are  shown  in  Fig.  3.  Table  III  shows 
that  the  estimated  values  corresponding  to  the  true  neighbor 
sets  are  close  to  the  true  values.  Note  that  when  extra 
neighbors  are  added,  the  corresponding  parameter  values  are 
very  small.  The  decision  statistics  corresponding  to  the 
models  considered  are  given  in  Table  IV  and  the  decision  rule 

(3.2)  picks  up  the  true  model. 

Fig.  3  shows  the  images  corresponding  to  the  different 
neighbor  sets  in  Table  III  using  an  identical  array  of  noise 
variables.  Windows  (1,2)  and  (1,3)  correspond  to  CMRF  models 
that  are  subsets  of  the  true  model  and  are  not  as  good  as  the 
original  in  (1,1).  Window  (1,4)  is  generated  by  the  true 
neighbor  set  with  estimated  parameters  and  is  very  similar  to 
(1,1).  The  pattern  (2,1)  corresponds  to  N^  =  {  (1 , 0) ,  (0 , 1)  } , 
which  has  no  common  neighbors  with  the  true  neighbor  set  and 
is  distinctly  different  from  the  original.  The  windows  (2,2), 

(2.3)  and  (2,4)  look  close  to  the  originals  since  the  generating 
models  include  the  original  neighbor  set.  But  the  decision 
rule  suggested  here  correctly  eliminates  these  models. 


Experiment  4:  Unilateral  approximations  to  CMRF  models 

Unilateral  approximations  to  CMRF  models  with  general 
neighbor  sets  of  dependence  are  of  interest  for  several 
reasons.  First,  the  estimation  of  parameters  in  unilateral 
spatial  autoregressive  models  can  be  handled  similarly  to 
one-dimensional  autoregressive  time  series  models.  Secondly, 
unilateral  neighbor  sets  make  possible  a  state  space  representa¬ 
tion  of  the  images,  which  is  useful  in  Kalman  filtering  of 
images  [13].  Though  it  would  be  futile  to  expect  that  "good" 
unilateral  approximations  can  be  found  for  any  arbitrary  CMRF 
model,  it  would  be  worthwhile  to  explore  the  adequacy  of  such 
approximations  even  for  some  specific  CMRF  models.  We  consider 
the  simplest  CMRF  model  with  =  {  ( 0 , 1)  ,  ( 1 , 0)  } .  As  observed 
in  [8]  several  unilateral  models  driven  by  white  noise  can  be 
constructed  in  increasing  order  of  complexity.  The  simplest 
unilateral  model  involves  dependence  on  the  north  and  west 
neighbors,  the  next  approximation  involves  dependence  on  the 
north,  west  and  southwest  neighbors,  and  so  on.  In  Fig.  4, 
we  have  shown  patterns  corresponding  to  several  such  approximate 
models.  Window  (1,1)  corresponds  to  the  original  pattern,  and 
(1,2)  to  the  same  neighbor  set  as  in  (1,1)  but  with  estimated 
parameters.  The  pattern  in  (1,3)  corresponds  to  the  unilateral 
neighbor  set  N  =  {  (-1 , 0) , (0 , -1)  } ;  the  parameters  are  estimated 
using  the  least  squares  method  and  the  pattern  was  reconstructed 
using  the  method  in  [5],  This  pattern  looks  surprisingly  similar 


to  (1,1)  and  (1,2).  A  careful  analysis  indicates  that  this 
is  to  be  expected.  This  unilateral  process  is  defined  such 
that  the  probability  distribution  of  y(i,j)  given  all  the 
observations  in  the  unilateral  region  is  equal  to  the  prob¬ 
ability  distribution  of  y(i,j)  given  y (i-1 , j ) ,y (i , j— 1) . 
However,  for  the  same  unilateral  neighbor  set,  the  conditional 
distribution  of  y(i,j)  given  all  others  except  y(i,j)  depends 
[8,21]  on  y(i-l, j) ,y(i+l, j) ,y(i, j-1) ,y(i, j+1) ,y(i-l, j+1) , 
and  y(i+l,j-l),  a  CMRF  model  that  "includes"  the  original 
CMRF  model  with  =  { ( 0 , 1) , (1 , 0) } .  As  illustrated 
in  Experiment  3,  the  patterns  corresponding  to  a  CMRF  model 
with  a  neighbor  set  N^and  another  CMRF  model  that  includes 
and  some  extra  neighbors,  look  very  similar.  Hence,  the 
simplest  unilateral  neighbor  set  N  =  {  (—1,0) ,  (0,-1) }  seems 
to  be  an  excellent  approximation  to  a  first  order  CMRF  model. 
Several  other  unilateral  approximate  patterns  and  bilateral 
patterns  corresponding  to  N  =  {  (0 , 1)  ,  (1 , 0)  ,  (0 ,  —  1 )  ,  (-1 , 0)  } ,  a 
SAR  model,  are  shown  in  Fig.  4.  The  details  of  the  models 
may  be  found  in  Table  V.  However,  in  more  general  situations, 
the  procedure  suggested  in  [1],  based  on  approximating  the 
spectral,  density  of  a  CMRF  model  by  successive  approximations 
of  unilateral  neighbor  sets,  should  be  used.  The  goodness 
of  approximation  can  be  visually  evaluated  as  described  above. 

By  comparison,  the  unilateral  approximations  with 
N  =  {  (0,-1),  (-1,0)  }  or  N  =  {  (0,-1),  (-1,0)  ,  (-1,1)  }  may  not  be 
good  for  SAR  models  with  N  =  {  (-1 , 0) , (1 , 0) ,  (0, 1) ,  (0 , -1) } , 


5.  Conclusions 


We  have  considered  some  aspects  of  statistical  inference 
methods  applied  to  two-dimensional  discrete  Markov  random  fields. 
Specifically,  we  have  considered  estimation  schemes  for  the 
parameters  of  the  Markov  model  and  decision  rules  for  the 
choice  of  an  appropriate  model.  We  have  illustrated  the  use¬ 
fulness  of  the  methods  by  using  synthetic  patterns.  Currently, 
work  is  in  progress  on  testing  the  estimation  schemes  using 


real  textures. 
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Table  I.  Details  of  CMRF  models  corresponding  to  Fig.  1. 
For  all  the  models  <*=30.00,  v -1.1111. 


Image 
(Fig.  1) 

Neighbor  set  N.^ 

Parameters 

(1,1) 

(1,0) , (0,1) 

ei,o=-2794' 

Yi=*1825 

(1,2) 

(-1,1) , (1,1) 

6-i,r* 28' 

ei,i=-14 

(1,3) 

(-1,1) , (1,1) 

e-i,i=-‘14' 

9.  .=.28 

J-  t  J- 

(1,4) 

(1,0) , (1,-1) , (0,1) , (1,1) 

Y<r-3357' 

9i,-r--25 

eo,i=-3246' 

0  =-.2126 

1  y  i- 

(2,1) 

(1,0) , (1,-1) , (0,1) , 

ei,o=-206^ 

ex  _1=.0536, 

(1,1) , (2,0) ,  (0,2) 

Yr*4467' 

6  x= .  0536  , 

e2,0=“'°123. 

0O,2=-0580 

(2,2) 

(1,-1)  , (1,1) ,  (2,0)  , 

Vr-23*!. 

0  =.4682 

J-  9 

(0,2)  ,  (2,-2)  ,  (2,2) 

e2,0=-°655' 

0O,2='°655 

®2,  -2~  "*°163, 

02,2="-°655 

(2,3) 

(-1,1) , (1,1) 

0-l,l=*28' 

01,1=-.22 

(2,4) 

(1,0)  ,  (0,1) , (2,0) , (0,2)  , 

9i,o=-12'  eo. 

1=~ ' 3-8  '  92,  0~ 

(3,0) , (0,3)  , (4,0)  , (0,4) 

60,2=-09' 

03 , 0=~ ’ 11 ' 

e0,3=-U' 

»,0— 0?,  °0,4 

(3,1) 

(1,0)  , (0,1)  , (3,0) ,  (0,3) 

0i,o=-16' 

eo,i=*10' 

03,O=-12' 

0O,3=-14 

(3,2) 

(1,0) , (0,1) , (3,0) ,  (0,3) 

9i,o='10' 

9o,r-16' 

03,O=-14' 

90 , 3=  * 12 

(3,3) 

(1,0) , (0,1) , (1,1) 

9i,o=-12' 

9o,r-10' 

(-1,1) , (3,0) , (0,3) 

0i,r-08' 

8_l,l=-.°9, 

03 ,  o~~ ' 11 ' 

90 , 3= ' 11 

I 


Parameters 


(3,4) 

(0,2)  ,  (2,0) 

90,2=*2794 ' 

02,OS‘1825 

(4,1) 

(1,0) , (0,1)  ,  (2,0)  ,(0,2) 

8i,o=,;L2' 

e0,r"-24' 

0  2 , 0= ‘ 18 • 

80 , 2~” ' 18 

(4,2) 

(1,0) , (0,1)  ,  (2,0)  ,  (0,2) 

61,0=~‘24' 

Vr*12, 

02,O='‘18, 

80 , 2~  ‘ 28 

(4,3) 

(1,0)  ,  (0,1) , (1,1) , (2,0) 

ei,o=~‘24. 

Vi-12 

(0,2) 

61,1=,U»  92 

,o=--18^  eo. 

(4,4) 

(0,1) ,  (1,0) 

9o,i=“-12' 

01,O=-18 

V 


Table  II.  Values  of  parameters  used  to  illustrate  the 
role  played  by  the  coefficients  in  synthesizing 
patterns,  u-30.00,  v*l. 11.11,  N^  =  f  (-1,1)  ,  (1,1) } 


Image  Parameters 

1Fiq-  21 _ ^ 


(1,1) 

.28 

-.14 

(1,2) 

.24 

-.14 

(1,3) 

.20 

-.14 

(1,4) 

.16 

-.14 

(2,1) 

.28 

-.10 

(2,2) 

.28 

-.06 

(2,3) 

.28 

-.18 

(2,4) 

.28 

-.22 

(3,1) 

.26 

-.20 

(3,2) 

.24 

-.18 

(3,3) 

.32 

-.14 

(3,4) 

.34 

-.16 

Table  III.  Details  of  models  fitted  to  the  synthetic 
data  generated  by  the  model  with  6_^  ^=-.14, 

0^  ^=.28,  a=30.0,  v =1.1111.  The  estimate  & 
is ' =30 . 005 . 


Image 
(Fig.  3) 

Neighbor  set  N^ 

0 

Estimates  of 
Coefficients 

(1,2) 

(1,1) 

1.1638 

6.  .=.3116 

(1,3) 

(-1,1) 

1.3520 

0_1  x=- . 2093 

(1,4) 

(-1,1) , (1,1) 

1.1033 

0  .  .=-.1410 

“X  t  X 

6^^  1=.  27875 

(2,1) 

(0,1) , (1,0) 

1.4934 

Vr--°o52 

0^^  q=-  .  0020 

(2,2) 

(-1,1) , (1,1) , (1,0) 

1.1033 

0_1  1=-. 14101,  0^^  .  27877 

ei,o=-'0051 

(2,3) 

(-1,1) , (1,1) , (0,1) 

1.1033 

0_x  1=- .  1410 ,  0X  ^  .27873 
0Q  L=-. 001717 

(2,4) 

(-1,1) , (1,1) , 

(1,0) , (0,1) 

1.1033 

0_1  1=-. 14101, 
x=. 27876, 

0i,o=--0049' 

0o,i=-.°009 

Table  IV.  Test  statistics  corresponding  to  the  fitted 
models  in  Table  III. 


Image 
(Fig.  3) 

Test  Statistic 

(1,2) 

1583.3 

(1,3) 

1637.2 

(1,4) 

1480.6 

(2,1) 

1676.3 

(2,2) 

1497.6 

(2,3) 

1497.1 

(2,4) 

1514.2 

Table  V.  Unilateral  approximations  to  the  first  order 

Markov  model  0,  =.2794,  0  =.1825,  a=30.Q,  v=l.llll. 

1  j  U  U  /  X 

The  estimate  3.  is  30.016. 


Image 
(Fig.  4) 

Neighbor  set 

0_ 

Estimates  of 
Coefficients 

(1,2) 

(-1,0)  ,  (0,-1) 

1.3117 

e-i,o='3634  '  Vr-2550 

(1,3) 

(-1,0)  ,  (0,-1)  ,  (1,-1) 

1.2991 

0_1  q= . 3575 ,  eQ  _1=.2166 

01  _i=*0944 

(1,4) 

(-1,0)  ,  (-1,-1) 

1.2991 

0  _1  Q= .  3575  ,  0_1  _1=.0004, 

(0,-1) , (1,-1) 

V-r-2165'  ei,-i=-0944 

(2,1) 

(-1,0)  ,  (-2,0)  ,  (-1,-1) 

1.2979 

6-l,0=‘3460'  0-2 , 0=  *  0291 ' 

(0,-1) ,  (-2,-1),  (1,-1) 

0_1  ^-.0033  ,  0Q  _1=.2159, 
0_2  _  .  0025  ,  01^_1=.O934 

(2,2) 

(-1,0), (-2,0),  (-1,-1)  , 

1.2976 

6-1, 0=‘ 3457  '  0-2,  d=-0285' 

(-2,-1),  (0,-1),  (0,-2), 

0_!  _1=-.0040,  0_2  ^=.0019, 

(1,-1) , (1,-2) 

e0,-f  *2140'  0O,-2=-0030 

ei,-l='0895,  91  ~2~ '  0127 

(2,3) 

(-1,0)  ,  (-2,0)  ,  (-1,-1), 

1.2974 

0_lfO=.3457,  0_2f(f.O287, 

(-2,-1),  (-1,-2),  (-2  , -2) , 

e_1  _]=-.0012,  0_2  _1=-.0038, 

(0,-1) , (0,-2),  (1,-1) 

0_1  _2=-.0120  e_2  _2=-.0028, 

(1,-2) 

0Q  _x  = . 2139 ,  0Q  _2= . 007 6 , 

61  _i=-0895'  ei  _2=-0132 

(2,4) 

(-1,0) , (1,0) , 

1.1582 

9-i,o=ei,o=-190°, 

(0,1) , (0,-1) 

0O,-l=eO,l=-1387 

(SAR  model) 

i 


Figure  1.  Examples  of  synthetic  generation 
of  patterns  using  CMRF  models  (for 
details  of  the  models  see  Table  I) . 


Figure  2.  Patterns  produced  by  models  with  the 
»'  same  neighbor  set  N,={ (-1,1) , (1,1) } , 

,  cx=3Q  .  0  ,  p=l.llll  but  with  different  sets 

X  of  values  for  the  coefficients  (see  Table  II) . 


Figure  3.  Synthetic  patterns  corresponding  to 

original  and  different  fitted  CMRF  models 
(see  Table  III  for  details  of  the  fitted 
models) . 


Figure  4.  Several  unilateral  approximations  to 
a  first  order  CMRF  model  (see  Table  V 
for  details  of  the  unilateral  models) . 
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